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Abstract: Neuronal brain activity in response to repeated stimuli can be per¬ 
ceived using functional magnetic resonance imaging (fMRI). In this paper, we 
develop a statistical model for fMRI data that estimates both the associated 
haemodynamic response function and the within and between trial variability 
through maximum likelihood estimation. We discuss our results in the context 
of other model-driven approaches, extending models already popular in the lit¬ 
erature, while removing the need for some of their assumptions. We consider an 
application to the motor cortex activity caused by a subject pressing a button 
and observe that the response changes signihcantly with task and through time. 
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1. Introduction 

In response to a stimuli, the activation of a brain region can be studied using functional 
magnetic resonance imaging (fMRI). Neural activity is detected by changes in the blood 
oxygenation level dependent (BOLD) haemodynamic response signal which is measured 
through monitoring magnetic resonance (MR) images of the subject’s brain taken at regular 
intervals. Psychologists and neuroscientists use fMRI experiments to identify MR signal 
intensity changes which correlate with the experimental paradigm presented to the subject. 
Inference can then be made regarding the location and time-course of the underlying neural 
activity. 

A typical fMRI experimental design consists of a number of trials, known as epochs. During 
each epoch, the subject receives a stimulus or performs a task corresponding to one of the 
experimental conditions. Each experimental condition is normally repeated several times 
and, for a large number of conditions, the experiment may be divided into a number of 
separate sessions. Repeating the experiment on several subjects enables inferences to be 
made about the population. 

At regular intervals during each epoch an MR scan gives a three-dimensional array of MR 
signal intensity measurements. Each entry in the array corresponds to a three-dimensional 
pixel, known as a voxel, in the image. The sequence of images taken during the experiment 
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leads to a four-dimensional data set. Physical constraints enforce a compromise between 
spatial and temporal resolution but modern scanners can typically record voxels of volume 
10 mm^ at intervals in the order of seconds. A complete brain image consists of several 
hundred thousand voxels. 

The aim of this paper is to build a statistical model for the haemodynamic response function 
(HRF) caused by a sequence of stimuli over a period of time and develop new inferential 
methods for fMRI data. Brain activation is highly variable, and establishing how the 
brain’s response changes to each stimulus (single-trial variability) is an important step to 
understanding the wider problem of how the different areas of the brain divide tasks and 
interact. 

The primary motivation for this work is the analysis of images from an ultra-high held MR 
scanner. The images are gradient echo planar images (EPI) taken at the Sir Peter Mansheld 
Magnetic Resonance Centre (SPMMRC), University of Nottingham, on the Philips 7T 
scanner. Scanners with an ultra-high held of 7T have only become recently available. The 
higher held increases the BOLD signal change and the signal to noise ratio, however, it 
also raises the confounding physiological noise levels. 

The plan of the paper is as follows. In Section 2 we introduce our case study and the exper¬ 
imental data acquired and examine previous work on this topic. In Section 3 we develop 
our statistical model for fMRI data and maximise the likelihood using an Expectation- 
Maximisation (EM) algorithm. In Section 4 we propose hypothesis tests for statistical 
inference on the parameter estimates. We present results for the case study in Section 5. 
Finally, we conclude with a brief discussion. 

2. Application and previous work 

2.1 Application 

Many fMRI studies involve a physical response to a visual or audio stimulus, for example 
selecting an option from a keypad. It is, therefore, vital for a wide range of applications 
that the time course of the BOLD signal change associated with neuronal activity, the 
haemodynamic response function (HRF), in the motor cortex can be modelled. For this 
reason, we acquired images of a volunteer’s motor cortex, which is an area of grey matter 
towards the posterior of the brain, that is responsible for controlling muscle movements. 

A dynamic scan was taken at 2 second intervals and each dynamic scan was composed of 12 
slices, each 3mm thick with a 0.7mm slice gap, taken sequentially at 1/6 second intervals. 
Each slice is composed of 64 x 64 pixels of size 1mm x 1mm. A total of ten visual stimuli 
were presented at 28.25 second intervals, which cued the volunteer to press a button once 
on odd-numbered trials or hve times on even-numbered trials. 

Acquired MR images are usually subjected to three preprocessing steps, namely slice¬ 
timing, realignment and smoothing, with the aim of removing confounding effects. Slice¬ 
timing corrects for the time difference between each slice being recorded within each dy¬ 
namic scan and transforms all the data within each scan to a single time point. Realignment 
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Figure 1: Schematic diagrams showing the difference between slice timing (top) and 
slice/trial timing (bottom) for a simplified paradigm with a dynamic scan of 3 slices and 
a 5.5 second interval between stimuli. Slice acquisition times are shown with a solid line, 
and the corrected timings with a dotted line. The stimuli were given at times marked “S”. 


corrects for the small movement in the volunteer’s position between each dynamic scan by 
translating and rotating each scan, in relation to the first, with a six-parameter rigid-body 
transformation. Finally, the images are spatially smoothed to improve the signal-to-noise 
ratio and allow for any remaining imperfections in anatomical alignment, at the expense 
of spatial resolution. The size of the spatial smoothing kernel should match the size of any 
potential activation we wish to detect. Typically, this is 1.5 times the voxel size, and our 
images were smoothed with an isotropic Gaussian kernel of 1.5mm full-width half-maximum 
(FWHM). 

We pre-processed our data using the SPM2 software package (Friston et al., 2002) but 
with one slight modification to the slice-timing. The slice-timing algorithm models the 
time series from each slice as a linear combination of sinusoids of different phases and 
frequencies and the data is then shifted forwards or backwards in time by effectively adding 
a constant to the phase of every frequency. Conventionally, each slice is shifted so that the 
time-series has the values that would have been obtained had the slice been acquired at 
the same moment as a reference slice. However, due to the “jitter” asynchrony between 
the stimulus presentation times and the scan acquisition rate, we implement slice-timing 
with a different time shift for each slice/trial combination so that the slices are corrected 
to the same post-stimulus time points for each trial. A schematic diagram of this can be 
seen in Figure [H 

During image preprocessing we automatically masked voxels at locations exterior to the 
brain such that they are not included in the following analysis. 
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2.2 Background 

Current statistical methods for fMRI data analysis tend to focus on localising areas of the 
brain where the response changes between two experimental conditions. The output is a 
classihcation in which each voxel is deemed active or inactive in processing the difference 
between the stimuli or tasks. A careful choice of experimental conditions will enable a 
localisation of the neurons involved in processing precise brain functions. There are two 
limitations to purely focusing on localising brain activity (Genovese, 2000). Firstly, fMRI 
data is correlated in space and time due to the underlying biological process and the 
measurement technique. For example, we would expect some correlation between voxels 
containing neurons performing similar tasks. Similarly the haemodynamic response will be 
influenced by a voxel’s proximity to major blood vessels, the location of which is unknown. 
Secondly, localisation does not inform the experimenter about changes in the response 
between trials or the temporal dynamics of the process. 

A variety of statistical techniques have been proposed for fMRI data. These can be divided 
into data-driven and model-driven methods (Lu et ah, 2005). Data-driven methods include 
independent component analysis (e.g. Beckmann and Smith, 2004), principal components 
analysis (e.g. Hansen et al., 1999) and cluster analysis (e.g. Goutte et ah, 1999). The 
advantage of these techniques is the ability to identify differences in the HRF between 
conditions without specifying a hypothesis or a model. Model-driven methods include the 
general linear model (e.g. Friston et ah, 1995b) and the deconvolution model (Goutte et 
al., 2000). Model-driven methods have gained in popularity due to the ease of interpre¬ 
tation and application. For example, the general linear model (GLM) approach can be 
implemented using the software packages FSL (Beckmann et ah, 2003), fMRIstat (Worsley 
et ah, 2002) and SPM (Friston et al., 2002). 

The majority of model-driven methods treat voxels independently, primarily for speed of 
computation. In the GLM method, for example, a linear model is htted with the voxel time 
series as the response variable and convolutions of a HRF with the experimental paradigm 
as covariates. The corresponding parameters, one for each condition, and contrasts of 
parameters can then be estimated at each voxel and plotted over a brain image to form 
a statistical parametric map (SPM). Hypothesis tests, such as t-tests and ANOVAs, with 
a null hypothesis of no activity can then be performed to localise voxels with signihcant 
activity. 

A basic implementation of the GLM method assumes the HRF shape is known. The 
standard canonical HRF employed by SPM2 (Worsley and Friston, 1995), for example, 
is shown in Figure [21 Its shape is derived from prior empirical evidence, demonstrating 
a positive peak 4-8 seconds after the onset of activity and a negative undershoot at the 
return to baseline, and is formed by the discretised summation of two curves taken from a 
gamma probability density function. The disadvantage of assuming a hxed shape is that 
the statistical tests have less power if the shape is incorrect. 

One solution to this problem is to increase the number of basis functions in the linear 
model. Studies have tried expanding the number of basis functions based on derivatives 
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SPM’s standard haemodynamic response function 



Figure 2: SPM2’s standard haemodynamic response function evaluated at two second 
intervals. 

with respect to time and dispersion (Friston et ah, 1998), principal component analysis 
of a large number of HRFs (Woolrich et ah, 2004b), terms of a Fourier series (Josephs et 
ah, 1997) and cosine functions (Zarahn, 2002). Taking the number of basis functions to 
the extreme, we could estimate the HRF shape by including one free parameter for every 
sampled time point, as implemented in the hnite impulse response (FIR) model (Glover, 

1999) . The disadvantage of these latter methods is that there are lower degrees of freedom 
and less power if the covariates are not orthogonal. 

A number of studies have tried more flexible models to estimate the HRF shape from the 
data. For example, the HRF shape has been parameterised as the discretised summation 
of inverse logit functions (Lindquist and Wager, 2007) or through cubic splines (Genovese, 

2000) . The disadvantage of these methods is that the parameters can usually only be esti¬ 
mated through non-linear algorithms or simulation techniques, such as Markov chain Monte 
Garlo (MGMG), which may have convergence problems and be computationally expensive. 
Other approaches include the selective averaging of responses to each experimental condi¬ 
tion (Dale and Buckner, 1997). From these estimates of HRF shape, inferences can then be 
made regarding the onset, peak latency and response duration. Both the GLM models and 
the more flexible models assume that some of the HRF parameters differ between voxels, 
although sometimes only the magnitude parameter, but not between trials. 

The assumption of independent trials is usually made because the BOLD signal change is of 
the order of 2.5% and the signal-to-noise ratio (SNR) is low. Averaging over trials, assuming 
that the response to each stimulus is dependent only on location and time since the last 
stimulus, improves the SNR. In truth, this is probably a false assumption. Experiments 
have shown that the HRF varies between voxels (Aguirre et ah, 1998) and between epochs 
(Duann et ah, 2002). Behavioural studies have shown “warm up” (Eysenck and Frith, 1977) 
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and “carry over” (Ward et al., 2001) effects, where peak performance is not achieved when 
starting a new task from rest or when switching immediately from one task to another. If 
a snbject makes a mistake in one trial, there will be a systematic effect on the behavionr 
in the next trial (Debener et al., 2005). Dnann et al. (2002) revealed “dramatic and 
nnforeseen haemodynamic response variations not apparent to researchers analysing their 
data with event-related response averaging and hxed haemodynamic response templates”. 

Adapting fMRI data analysis techniqnes to model trial-to-trial variability within a voxel 
is non-trivial. In the GLM context, it can be done nsing a non-additive analysis of vari¬ 
ance (ANOVA), where the time since the stimnlns and the epoch nnmber are taken as 
the treatment gronp and block effect respectively (Anffermann et ah, 2001). An epoch- 
specihc parameter scales the treatment effect, corresponding to a change in amplitnde over 
trials, and the parameters are estimated using maximum likelihood estimation (MLE). Al¬ 
ternatively, an extension of the FIR model estimates a separate HRF function for each 
experimental condition, which is allowed to vary in amplitude over trials within the con¬ 
dition (Hinrichs et ah, 2000). The HRFs and the amplitudes are estimated through a 
non-linear optimisation algorithm assuming a Toeplitz-type covariance structure that is 
constant over voxels. 

Devising a model which allows more complicated variation between trials greatly increases 
the number of parameters. For example, an alternative extension of the FIR model is to 
estimate an HRF for each trial but constrain the estimate to lie in a neighbourhood of the 
original HRF estimate using non-linear optimisation techniques (Lu et ah, 2005). Another 
approach is to model the response to each trial using splines with parameters for lags and 
amplitudes (Genovese, 2000). This model is very flexible and, if prior distributions are 
specihed for the parameters, the posterior distribution can be maximised or sampled from 
using MGMG. 

Our approach will be to use a linear mixed effects model with an overall HRF, with one 
free parameter for each time point, a single hxed effect scale parameters for each voxel 
and, hnally, random effects for errors within and between epochs. Our approach provides 
a hexible model with a manageable number of parameters. 

3. Modelling the data 

Let Y[t) be the observed MR signal from a specihc voxel at time t = 1,... ,N, where 
N is the total number of images in the experiment. The MR signal intensity is affected 
by covariates such as movement, blood how and oxygen levels. Therefore, let X be the 
corresponding Nxq design matrix, composed of image preprocessing and physiological data, 
with the mean of each column being zero. The response variable and the design matrix 
are both subjected to temporal hltering using the hrst few non-constant basis functions 
of a one-dimensional cosine transform to remove low-frequency confounds, such as MR 
scanner drift. This is a standard SPM2 preprocessing step. The mean MR signal intensity 
is also subtracted from the data to set the voxel’s baseline signal to zero. Following this 
transformation, we denote the design matrix as X and let be the haemodynamic response 
of the Ah voxel. The pre-processed data for two voxels are shown in Figured 
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Pre-processed data for two voxels 



Figure 3: The pre-processed data from two example voxels. The stimuli were presented at 
the times marked with dotted lines. One shows signs of activation (top), whilst the other 
is mainly noise (bottom). 

Let E be the number of epochs (or trials) in the experiment and let T be the number of 
images during each epoch (i.e. N = T x E). Let h{t) be an unknown vector of length T 
that represents the haemodynamic response in each epoch. We expect the response at each 
stimuli to be similar, so we let fi = lE®h represent a convolution between h{t) and stimulus 
onset. Let Eg and Et denote variability between and within events, respectively, where the 
subscript denotes the dimension of the covariance matrix. We cannot estimate /r. Eg and 
Et’ for individual voxels due to the number of parameters involved, so we introduce some 
regularisation across voxels. Figure [3] shows that not all voxels within the brain are active 
during an experiment, so we base the estimates on active voxels only. This is similar to 
SPM2’s method for estimating variability based on voxels demonstrating significant regres¬ 
sion (Mumford and Nichols, 2006). We assume, however, no prior information regarding 
which voxels are active. Therefore, we formulate a multivariate Gaussian mixture model 
for two classes of activity (active and inactive). 

We model the response for voxel i, with 

active component: T) ~ Ei = E^; 0 E-r) and 

inactive component: Yi NTE{i ^2 — Xbi, E 2 — cr^ln)- 

where A denotes the scale of the response at the ith voxel, 6* is the parameter vector for the 
design matrix and cr^ is a scalar. To remove the arbitrary scaling of the response parameter 
we constrain the haemodynamic response function to have unit norm, i.e. ||h|| = 1. Note 
that the model for the active component is the mixed effects model, AjA: = A^fc + + 

Uijk, j = 1, • • •, -E, A; = 1,..., T, where 

(ujii,..., UietY' ~ N (0, Ee 0 E-r), 

A, h and bi are hxed effects and Uijk are the random effects. The Gaussian mixture model 
is then /(F*) = p/i(A) (1 -p)/ 2 (A), where /i ~ iV(/ii,Ei) and /2 ~ iV(/i 2 ,E 2 ). The 


7 
























log-likelihood for the data is, 


/(p,/il,/i 2 ,Si,S 2 |F) = X^log [p/l(yi|/il,Si) + (1 -p)/ 2 (Fi|/i 2 ,S 2 )], 

i=l 

where V is the number of voxels. 

Maximising this likelihood function is non-trivial. A solution does exist if we are prepared 
to constrain Si = S 2 (Day, 1969) but this requires the false assumption that both classes 
of activity demonstrate the same covariance structure. An alternative solution is to assume 
each voxel is either active or inactive and then evaluate the likelihood for all possible voxel 
allocations. This is equivalent to minimising |Si|^i |S 2|'^^5 where Sj is the MLE of the 
covariance matrix for the Vi voxels in group i = 1,2. However, there are possible 
combinations, which is computationally impractical. 

Instead, we utilise the Expectation-Maximisation (EM) algorithm (Dempster et ah, 1977) to 
maximise the likelihood and estimate the model parameters from the data. The algorithm 
requires that the data is augmented with latent variables, which are effectively missing 
data. Let Zi be a Bernoulli random variable with P(Zi = 1) = p. Let V fi if Z, = l 
and 17 ~ /2 if = 0. The joint density of {Zi,Yi) is 

/(Z„K.) = |p/,(K,y‘ 1(1 -p)/2(K,)]'-^' , 

and the log likelihood for the data is, 

/(p,/il,/i2,Si,S2|2', F) = 

V 

Y, Z^ log [p/i(F.|/ii, El)] + (1 - Z,) log [(1 - p)/ 2 (F,|p 2 , S 2 )]. 

i=l 


The EM algorithm is an iterative procedure with each iteration composed of expecta¬ 
tion and maximisation steps. The expectation step estimates the values of the latent 
variables and then the conditional log-likelihood is maximised by differentiation. Let 
0 ' = {p, /ii, p 2 , El, E 2 ) be the current parameter estimates and let 0 denote the parameter 
values at the next iteration then. 


Q(0|0',F) =E[/(0|Z,F)|F,0'], 


V 


Y^(z^^y,eVog 


2=1 


p/i(r.|6)) 

L(i-p)/2(y.ie) 


+ log 1(1 - p)/2(r.|e)], 


where Q has its usual EM interpretation and, 
E(Z,jV,0') = f(Z, = IjV, O') = 


p/i(D.| 0 ) 


p/i(F,|0) + (1-p)/2(F.|0) 


= Pi- 



Due to the small values of fi and / 2 , it is computationally more feasible to calculate 
Pi = 1/(1 + e^), where c = log(l - p) - log(p) + log(/ 2 ) - log(/i). 

We write the function, Q, in terms of our model parameters, 

V 

Q(0|0',F) = ^pilog[p/i(Fi|/ii,Si)] + (1-pi)log[(l-p)/ 2 (l/|/i 2 ,S 2 )], 

i=l 

VET f , / \ PiE , I „ I PiT 

= -^log7r + ^|pilog(p) - — loglSrl - — loglSgl 

i=l 

-|(y, - /iA - xkfE^^Yi - fip, - xh) 

+ (1 - Pi) log(l - p) - log(a2) 

- Xhi)'^{Y, - Xhi)y 


The model parameters can be estimated in turn by iteratively maximising Q conditional 
on the current values of the other parameters. The conditional MLEs are given below, with 
their derivation provided in the appendix. 


p 


A 

h 

h 




a 


V 


vY.P‘’ 


i=l 

{A^A)-^A^ - M) + (1 - , 

E* Ei Y.kPYjkPi(Xij - Xjbi 

{Y^iPif^i) (EjEfcejfc 

1 V' 

^^PiRiTiE Ri , 


EYtlP^ 

. V 


Pitt 

1 

ET^,(l-p.)^ 


^PiR^E^^Ri, 


V 


^(1 - p,)(Y, - Xh)^{Y, - Xh). 


where A = pjX^S/^X + (l—pj)X^S^^X, Yij and Xj are the rows of 1/ and X, respectively, 
corresponding to event j, Cjk is the j, kth entry of and Ri is the TxE matrix of residuals 
in the active model. 


EM algorithms are sensitive to the starting estimates of the model parameters. For our 
application we initially set h to be the standard haemodynamic response function shown 
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in Figure m and set S® and to be identity matrices. A reduced version of the model 
that assumes all voxels are active was then htted. The model parameters for the reduced 
model can be estimated using a simpler version of the algorithm presented above but with 
no expectation step, t-tests of the resulting f3i parameters provided an initial classihcation 
of active and inactive voxels for the EM algorithm. This classihcation was used to produce 
starting estimates of the proportion of active voxels, p, and the covariance matrices for the 
two groups, El and E 2 . 

4. Statistical Inference 

Our choice of model suggests that active voxels will have pi ~ 1 and a large positive 
f5i parameter. Active voxels, therefore, can be identihed by testing the null hypothesis, 
Hq : /3i = 0, versus the alternative hypothesis. Hi : (3i > 0, for i = under the 

model, Yi ~ + -A6j,Ei). A transformation gives Y* ~ NTE{f^*/3i + H*bi, Ite) 

where Y* = p* = and X* = The test statistic under TE — q — 1 

degrees of freedom is, 

^ ~ ^TE-g-1, (1) 

where Sf = ||F/ - p*A - X*bif/{TE -q-l). 

There has been much debate in the neuroimaging community regarding how to set the 
statistical signihcance level in fMRI studies (Marchini and Ripley, 2000), due to problems 
in estimating the magnitude of the response and the correlation in spatially and temporally 
correlated data. SPM2 utilises set-level inference, using distributional approximations from 
the theory of Gaussian random helds (Friston et ah, 1996). The method assesses the 
probability of obtaining c or more clusters, containing v or more voxels. Our voxel-by- 
voxel approach is the simplest case of this, allowing clusters of just one voxel. The use of 
Gaussian random helds has potential advantages over SPM2’s previous methods (Friston et 
ah, 1994; Friston et al., 1995a; Worsley and Friston, 1995), which corrected the signihcance 
level for temporal correlations only. Our method of pre-whitening the data, using the 
estimated covariance matrix to transform from Y, p and X io Y*, p* and X*, is similar 
to that of Worsley et al. (2002), although they restrict the covariance matrix to be of 
auto-regressive form. We correct for multiple comparisons by adjusting the initial p-value 
threshold of 0.001 to control the false discovery rate (Benjamini and Hochberg, 1995) using 
the adaptive method (Benjamini and Hochberg, 2000). 

The goal of single-trial variability analysis is to identify trends in the way active voxels 
respond in space and time to stimuli. Here a key ingredient of our work is to examine the 
random ehects at the active voxels. We carry out this task by investigating the principal 
components scores of the htted random ehects covariance matrices at the active voxels. 
Principal components (PG) analysis of will highlight key diherences in the response of 
voxels to diherent events. The fcth PG score for voxel i at event j is, 

^ijk Tfc 
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where 7 ^ is the fcth PC loading. We model the kth PC score using an analysis of variance 
with event number as a categorical variable. 

5. Results 

In this section we apply statistical models to the motor-cortex data in Section 2. The post¬ 
stimulus times chosen for our study were |, 2|, 4|,..., 26| seconds, as these correspond to 
the middle-slice acquisition times during the hrst trial, so T = 14. There were E = 10 
epochs, so N = TxE = 140. We choose q = 6 , using the six realignment parameter vectors 
as covariates because physiological data were unavailable. Recall that X is the N xq design 
matrix. The model proposed in Section 3 has a much more complicated structure than the 
simpler linear model currently implemented by SPM2, so we compare 5 models of varying 
complexity. The hrst two linear models we consider are. 

Model 1: R* ~ NTEiPif^s + Xbi, Ite)'-, 

Model 2: Yi ~ NTEWifJ^ + Xbi, ct'^Ite)] 

where is the HRF used by SPM2 and p is an HRF estimated from the data. In both 
models the shape of the HRF is constant for each voxel and trial but the magnitude of the 
response varies at each voxel through the estimation of /5j, i = 1, ...,V. We compare these 
to Gaussian mixture models where the active component is modelled by. 

Model 3: Active component Yi ~ NTE{,(3iIi> + Xbi, ^e ^ It)] 

Model 4: Active component Yi ~ NTEi.f^il'' + Xbi, Ie® ^t); 

Model 5: Active component Yi ~ NTEif^ifJ^ + Xbi, ^e ® ^t)] 

and, in each case, the inactive component is modelled by Yi NrEiXbi, a'^IxE)- In these 

models the expected response and covariance matrix structure is different for the active 
and inactive components, and the models only differ in the complexity of the covariance 
structure for the active component. 

For models 3-5, the EM algorithm was run until the change in the norm of parameter 
estimates between successive iterations fell below a specihed tolerance of 0.0001. The EM 
algorithm took approximately 50 iterations to converge for the model and data used in 
this study. The algorithm only takes a few minutes to compute but it will require more 
iterations if the amount of data or covariates is increased. 

The Akaike information criterion (AIC) and Bayesian information criterion (BIG) were 
calculated for each of the hve models using the log-likelihood, logL, the number of model 
parameters, P, and the number of observations, n. The results are shown in Table [H 
The more complicated models are favoured by a low AIG and BIG due to the addition 
of relatively few parameters. The reduction from model 1, based on SPM2 analysis, to 
model 2 highlights the beneht of allowing the data to estimate the haemodynamic response. 
Similarly, the large increase in the log-likelihood from model 2 to models 3-5 shows the need 
for separating the active and inactive components with a mixture model. A mixture model 
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Figure 4: A scatter plot of pi versus (3i (left) and a histogram of (3i values (right), showing 
active voxels (red) and voxels that share their covariance structure (blue). 


allows the variability of active responses to be estimated using only an appropriate subset 
of the voxels. The reduction in AIC and BIC between model 5 and model 4 suggests 
an improvement in the relative model fit and gives evidence that responses to sequential 
epochs are correlated and should not be treated as independent responses. Based on Table 
[H it seems reasonable to model the data using model 5. 


Table 1: The deviance and AIC for five models. 


Model 

P 

logL 

AIC 

BIC 

1 

70434 

-12481601 

25104069 

25940417 

2 

70447 

-12348551 

24837996 

25694502 

3 

80566 

-3898027 

7957185 

8936720 

4 

80616 

-3894722 

7950676 

8931818 

5 

80671 

-3889950 

7941242 

8922053 


A scatter plot of pi versus /9j for model 5 is shown in Figure 01 It highlights that voxels 
generally fell into three categories. Firstly, voxels with pi ^ 1 and Pi significantly greater 
than zero by the t-test in Equation ([1]), are the voxels which are truly active. Secondly, 
voxels with Pi ~ 1 and with /?, not significantly greater than zero are voxels which share a 
covariance structure with the active voxels but do not demonstrate a large haemodynamic 
response. Finally, voxels with pi ^ 0 have an MR signal that consists mostly of white 
noise. This provides some justification for using the Gaussian mixture model since truly 
active voxels nearly always have p* ~ 1 and we have removed some of the influence of voxels 
containing noise from the estimates of variability between and within trials. 

The voxels where Pi is significantly greater than zero are highlighted on the activation map 
in Figure |5l Note that head movement during scanning has caused part of the image in slice 
12 to not be consistently recorded throughout the experiment and consequently this portion 
of the image is masked during analysis. Slices 1-6 are not displayed for clarity and because 
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they contain less activation. Cluster analysis of the voxel co-ordinates revealed three large 
groups of voxels, colour coded in Figure O These areas of activation are similar to those 
detected using SPM2 with the simplest model, see Section 6, but here they illustrate where 
single trial variability is assessed. 

Figure [7] shows that the first three principal components of Sy explain a litle more vari¬ 
ability than the others. The plots show the mean response, /dp, and the effect of the first 
three principal components of Sy on the mean, /3p ± where /? = 10 and A^, 'jk are 

the kth eigenvalue and eigenvector of Ht- Voxels with a low PCI score show a much earlier 
rise in the response and a slightly later decay. The spatial distribution of PCI scores for 
Slice 8 is shown in Figure [H PC2 accounts for the width of the response, with low scores 
corresponding to a longer period of activation. Voxels with a high PCS score have a much 
later peak and a stronger response. 

The PC scores were modelled using a two-way analysis of variance (ANOVA) with event 
number and cluster as categorical factors. The fitted values from the ANOVA model are 
shown in Figure O Both PCI and PCS display a strong association with event type. 
There are more similarities between cluster 1 and 2 responses compared to the cluster S 
response. This mirrors the spatial relationship of these clusters and could reflect differences 
in each region’s function in processing the stimuli or task. For example, both motor and 
somatosensory areas are involved in this task and are captured in the MR image. 

Changes in haemodynamic response for a particular cluster can also be seen through plot¬ 
ting the htted values. 


3 

^cj Pck' T ^ ^ ■^cjk'^ki 
k=l 

where (3c is the mean value of (3i for the voxels in cluster c and Scjk are the htted PC scores 
from the ANOVA analysis. Figure [10] shows the htted values which have been interpolated 
using a cubic smoothing spline to aid interpretation. The plots clearly show that for clusters 
1 and 2 the response is much stronger and peaks later with a larger undershoot in trials 
where the volunteer presses the button multiple times. This ehect is less pronounced in 
cluster 3. However, it should also be noted that the one-press trials generally start from 
a lower baseline caused by a prolonged undershoot in the preceding multiple-press trials. 
There is a tendency for the response to be largest in the middle of the experiment. 

A principal component analysis of Eg revealed that the loadings do not display much 
structure. However maximum likelihood estimation of shows a relatively large increase 
in the likelihood compared to E^ = Ie and so it is worthwhile including this term in the 
model. 

6. Discussion 

In this paper we have developed a statistical model for analysing single trial variability in 
fMRI data. The chosen model produced a signihcantly higher likelihood than other simpler 
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Slice 7 


Slice 8 



Figure 5: The locations of active voxels in the motor cortex. Lighter colours indicate 
increased significance. 
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Figure 6: The location of active clusters, projected into the x-y plane. 
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Figure 7: The hrst plot shows the percentage of variability in each principal component 
of St’. The remainder show 10 times the mean response plus (green)/minus (red) the hrst 
three principal components of Ht- 
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Figure 8: PCI scores for voxels in Slice 8 for the ten different events. Odd-numbered 
trials required single button-presses (top row) and even-numbered trials required multiple 
button-presses (bottom row). PC scores ranged from -10 (red), through zero (blue) to -|-10 
(green). 


models, see Table [H including those implemented by software packages currently available. 
For comparison, we analyse the same data using two current statistical software packages, 
MELODIC (Beckmann and Smith, 2004) which implements a probabilistic independent 
components analysis, and SPM2 (Friston et ah, 1995b), which hts a linear model at each 
voxel, with a design matrix composed of covariate information and a basis function of a 
typical haemodynamic response. These methods are data and model driven, respectively, 
and their resulting activation maps are shown in Figures [H] and [121 Unsurprisingly, the 
active areas for both these methods have similar location to the active voxels detected using 
our model. We can conclude, therefore, that our model indicates areas of the brain active 
in performing the task of pressing a button, which concur with other methods of analysis. 

The key advantage of our method, however, is that we can also obtain maximum likelihood 
estimates of the shape of the response and the within and between trial variability. The 
three main sources of variation, displayed through principal components analysis, were 
found to be the timing of the initial rise in response, the length of the response and the 
strength of the peak. The analysis of principal components scores through a second-level 
model is similar to the multi-level models used in analysing between subject variability. 
fMRI group analysis for a particular voxel combines a low-level hxed effects model, similar 
to model 1 above, for the within-subject analysis with a high-level random effects model 
for the between-group analysis (Woolrich et ah, 2004a). 

For example, the low-level model for the ith voxel in the kth subject is Yik = Xkfiik + Cifc, 
where Tjfc is the recorded MR signal, is the design matrix composed of the haemody¬ 
namic response, experimental conditions and physiological covariates and eik r'U jv(o.4u). 
The parameter vector is given by the generalised least squares solution after pre-whitening 
of the data. The correlation matrix 14 can be estimated but must be spatially regularised to 
avoid biased estimates. The low-level model, therefore, is similar to the active component 
of model 5 but we have, additionally, estimated the haemodynamic response. 
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Figure 9: The fitted values and standard errors from the ANOVA analysis of PC scores. 
The colours correspond to cluster 1 (red), cluster 2 (green) and cluster 3 (blue), as shown 
in Figure [HI The results for each cluster have been offset to give greater clarity. 
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Figure 10: The fitted haemodynamic response for the three clusters for one-press trials (left 
column) and hve-press trials (right column). Individual trials are represented by different 
shades, with darker colours indicating later trials. 
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Figure 11: Highlighted voxels indicate locations where the signal is present. 



Figure 12: Highlighted voxels have signihcant non-zero activity under SPM2’s model with 
a p-value below 0.001. 
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In the multi-level group study approach, the parameter vector, f3ik, for each subject is taken 
as the response variable in the between-subject model, [Pfi ■ ■ ■= Xg(3g + eg where 
Xg is the design matrix composed of subject level covariates such as gender, age, health 
condition etc. No frequentist solution to these equations exist if the variance components 
are unknown so the summary-statistics approach is to solve the low-level model for each 
subject hrst and then solve the high-level model based on these estimates (Mumford and 
Nichols, 2006). 

Theoretically, this model could be applied to single-trial variability analysis by estimating 
the parameter vector for each voxel. Pi, at each trial and then applying the high-level 
model to analyse the between-trial variability. In other words, consider more fixed effects. 
This approach may not yield robust estimates given the limited data and high number of 
parameters unless the high-level model was applied to different sessions of the experiment 
with the same subject. This, however, requires the same subject to participate in the 
experiment multiple times to investigate consistent within-session variability (Lion et ah, 
2006) which may not be practically, or ethically, possible. 

Our model, however, could easily be extended to cope with multiple sessions or subjects. 
At present, we analyse the PC scores by looking for trends across voxels within a cluster. 
This is not an optimal approach because neighbouring voxels are likely to be correlated 
but with more runs or subjects, we could invoke a voxel-by-voxel approach to the PC 
analysis. It may also, depending on the application, be appropriate to treat the epoch 
number as a continuous variable and model the PC scores with a more general linear 
model. This approach was not implemented here due to the limited number of epochs of 
each experimental condition and it would not reveal interesting non-linear changes. 

Our lower-level model could also be adapted for multiple sessions or subjects. One variation 
would be to incorporate a different haemodynamic response function, h, for each subject or 
experimental condition. We would suggest, however, that including more basis functions 
would generate a larger number of possible contrasts and make physiological interpretation 
harder. The current model captures the inter-trial variability in the covariance model alone. 

It remains unclear, however, if the observed single-trial variability is caused by changes 
in neuronal processes or by natural variation in the reaction time of the volunteer and 
the strength of the button press. The amount of noise in the data could be reduced with 
the inclusion of physiological data as covariates but these data were unavailable for the 
experimental results presented here. A further rehnement of this model might also include 
an underlying continuous space/time model, using the real time and location that each 
slice was recorded and, therefore, eliminating the need for image preprocessing. 

The experiment paradigm in our application, with single and multiple button presses, 
also raised interesting questions regarding the relationship between neuronal activity and 
haemodynamic response. It is well documented that this relationship is non-linear, but 
predicting the change in response to a longer period of activation, for example, is less 
well understood (Mechelli et ah, 2001). SPM2 adjusts its basis function by convolving the 
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period of neuronal activity with its standard haemodynamic response function, using the 
Volterra series (Friston et ah, 1998). Independently, the balloon/Windkessel dynamical 
model of how the haemodynamic response is influenced by the underlying physical changes 
in the blood vessels was developed (Buxton et ah, 1998). The model suggests that increased 
blood flow inflates a venous “balloon”, diluting and expelling deoxygenated blood causing 
an increase in the BOLD signal. As the flow decreases, the balloon deflates reducing 
the discharge and increasing the concentration of deoxygenated blood, causing the post¬ 
stimulus undershoot. These models have been shown to be consistent, in that the Volterra 
kernels which best represented those derived from empirical evidence, also had biologically 
plausible estimates for the balloon/Windkessel model parameters (Friston et ah, 2000). 
In this study, we have developed a model that distinguishes between different responses, 
caused by different neuronal activity, without making assumptions regarding the underlying 
non-linear relationship. 

The methods and results given in this paper obviously pertain to one experiment carried 
out with one volunteer. It would be interesting in further work to apply the model and 
methodology to other volunteers, or repeat it with the same volunteer, to examine if the 
sources of trial variability found in this study are common to all subjects or experiment 
repetitions. If this proves to be the case, it will greatly enhance the ability for neurolo¬ 
gists to reach conclusions on voxel activation where a traditional repeated trial experiment 
paradigm is either impossible to conduct or impractical with the resources available. 

To conclude, in this paper we have proposed a statistical model to examine single-trial 
variability for analysis of ultra-high held fMRI. This model is an extension of the SPM2 
model which is very popular in the literature. Importantly, through examining the PC 
scores, we are able to show how the response changes with the task and through time and 
demonstrate that these changes are statistically signihcant. The model also estimates the 
shape of the mean response and, for the case study considered, is similar to that of prior 
studies but displays a greater undershoot than previously reported, particularly when the 
peak is stronger and wider in the cases where the volunteer presses the button multiple 
times. 
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Appendix 

We estimate the model parameters in turn, maximising Q conditional on the current values 
of the other parameters by setting the derivative of Q with respect to the parameter equal 
to zero and solving for the parameter. Firstly, we consider the proportion of active voxels, 

p. 


dQ ^ Pi (1 - Pi) \ ^ ^ 

Qp ^\p (1 - p) j 



2 = 1 


For the scale parameter at each voxel, i = 1,..., U, 

|| = -P.(Ap^srV-h^sri(u-x6,)) = 0, 

^ {jxfii - {Yi - Xhi)) = 0. 
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Hence, 


A = - xh). 


For the parameter vector, 

^ {Xk - (A - /iA)) + (1 - p,)x^^^\xk - A) = 0. 

uhi 

Hence, 

k = {A^A)-^A^ {p.X^E^^Y, - /iA) + (1 - P^)X^^^^Yi) , 

where A = piX'^T.^^ + (1 — pijX'^T.^^■ Let Aj be the response from voxel i at event j, Xj 
the rows of X corresponding to event j, and let ejk be the entry at the jth row and kth 
column of then differentiating Q with respect to h gives, 

^ f V E E "I 

Piejk{Yij - hl3i - Xjbi)i:j,^{Yik - h/Si - Xkh) I 
^ t i=l i=l k=l j 

V E E 

E E Z (PV> - - xfiP) 

i=l j=l k=l 

Hence, 

^ (E* Ei Efc PYjk/^i(Yij - Xjbi)^ 

iJ2iPil3!) (EjEfcejfc) 

The relative sizes of A and p = 1 e ^ h are arbitrary so we can artihcially rescale them at 
each iteration such that ||h|| = 1 without changing the value of the likelihood. Let Ri be a 
T X E matrix of residuals, where the jth column is Aj — hA — Xjbi, and let A = Ri^E^Rf- 
Then, 


= 0 , 

= 0 . 


dQ 




1log I I -PAt^ 


= 0 . 


Given, (9log |Sy |/i9S.j, = 2Sr — Diag(ST), and (9tr(S^^S')/(9S^ = 2S — Diag(S'), if 


and S are symmetric, then, 

y 


5^teF;[2ST-Diag(ST)]-pA2^*-Diag(A)]} = 0, 
1 


i=l 


V 


pEj’ — Diag(S7’) 


EEtiPi ^ 


- Diag(A)] = 0 . 
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Let M = Tjrp — then we require 2M — Diag(M) = 0. Hence M = 0 and, 


V 


St — 


EEtiPi 


^PiSi 


By analogy, let T* = then. 


V 


St = 


^EtiPi lit 


^PiTi. 


Finally, 


Hence, 


V 


3Q ^ r-d - 


da^ ^ 1 2 ( 7 ^ ' 2 ( 7 ^ 

i=l 


= 0 . 


V 


= 


ETYEU^-P^)^, 
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